import math

a = 6378137
G_M = 398600441800000
omega = 0.00007292115
f = 1 / 298.257223563
b = 6356752.31414
e_prime_2 = 0.00673949674228
E = math.sqrt(a * a - b * b)
u = b
e_prime = E / b


m = omega * omega * a * a * b / G_M
U0 = (
    G_M * (1 - e_prime_2 / 3 + (e_prime_2 * e_prime_2) / 5) / b
    + omega * omega * a * a / 3
)

q0 = 0.5 * ((1 + 3 * b * b / (E * E)) * math.atan(E / b) - 3 * (b / E))
q0_prime = 3 * (1 + u * u / (E * E)) * (1 - (u / E) * math.atan(E / u)) - 1
gama_a = G_M * (1 - m - (m * e_prime * q0_prime) / (6 * q0)) / (a * b)
gama_b = G_M * (1 + m * e_prime * q0_prime / (3 * q0)) / (a * a)
phi1 = (23 + 12 / 60 + 23 / 3600) * math.pi / 180
phi2 = (-23 - 12 / 60 - 23 / 3600) * math.pi / 180
h1 = 2000
h2 = 6000

gama1 = (
    a * gama_a * math.cos(phi1) * math.cos(phi1)
    + b * gama_b * math.sin(phi1) * math.sin(phi1)
) / math.sqrt(
    a * a * math.cos(phi1) * math.cos(phi1) + b * b * math.sin(phi1) * math.sin(phi1)
)
gama2 = (
    a * gama_a * math.cos(phi2) * math.cos(phi2)
    + b * gama_b * math.sin(phi2) * math.sin(phi2)
) / math.sqrt(
    a * a * math.cos(phi2) * math.cos(phi2) + b * b * math.sin(phi2) * math.sin(phi2)
)
gama1_h1 = gama1 * (
    1
    - (2 / a) * (1 + f + m - 2 * f * math.sin(phi1) * math.sin(phi1)) * h1
    + 3 * (h1 * h1) / (a * a)
)
gama2_h1 = gama2 * (
    1
    - (2 / a) * (1 + f + m - 2 * f * math.sin(phi2) * math.sin(phi2)) * h1
    + 3 * (h1 * h1) / (a * a)
)
gama1_h2 = gama1 * (
    1
    - (2 / a) * (1 + f + m - 2 * f * math.sin(phi1) * math.sin(phi1)) * h2
    + 3 * (h2 * h2) / (a * a)
)
gama2_h2 = gama2 * (
    1
    - (2 / a) * (1 + f + m - 2 * f * math.sin(phi2) * math.sin(phi2)) * h2
    + 3 * (h2 * h2) / (a * a)
)


print(gama1)
print(gama2)
print(gama1_h1)
print(gama2_h1)
print(gama1_h2)
print(gama2_h2)
